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Protoplanetary gas disks are likely to experience gravitational instabilites (GI's) during 
some phase of their evolution. Density perturbations in an unstable disk grow on a dynamic 
time scale into spiral arms that produce efficient outward transfer of angular momentum and 
inward transfer of mass through gravitational torques. In a cool disk with rapid enough cooling, 
the spiral arms in an unstable disk form self-gravitating clumps. Whether gas giant protoplanets 
can form by such a disk instability process is the primary question addressed by this review. 
We discuss the wide range of calculations undertaken by ourselves and others using various 
numerical techniques, and we report preliminary results from a large multi-code collaboration. 
Additional topics include - triggering mechanisms for GI's, disk heating and cooling, orbital 
survival of dense clumps, interactions of solids with Gl-driven waves and shocks, and hybrid 
scenarios where GI's facilitate core accretion. The review ends with a discussion of how well 
disk instability and core accretion fare in meeting observational constraints. 



1. INTRODUCTION 

Gravitational instabilities (GI's) can occur in any region 
of a gas disk that becomes sufficiently cool or develops a 
high enough surface density. In the nonlinear regime, GI's 
can produce local and global spiral waves, self-gravitating 
turbulence, mass and angular momentum transport, and 
disk fragmentation into dense clumps and substructure. The 
particular emphasis of this review article is the possibility 
{Kuiper, 1951; Cameron, 1978), recently revived by Boss 
(1997, 1998a), that the dense clumps in a disk fragmented 
by GI's may become self-gravitating precursors to gas giant 
planets. This particular idea for gas giant planet formation 
has come to be known as the disk instability theory. We 
provide here a thorough review of the physics of GI's as 
currently understood through a wide variety of techniques 
and offer tutorials on key issues of physics and methodol- 
ogy. The authors assembled for this paper were deliberately 



chosen to represent the full range of views on the subject. 
Although we disagree about some aspects of GI's and about 
some interpretations of available results, we have labored 
hard to present a fair and balanced picture. Other recent re- 
views of this subject include Boss (2002c), Durisen et al. 
(2003), and Durisen (2006). 

2. PHYSICS OF GI's 
2.1 Linear Regime 

The parameter that determines whether GI's occur in thin 
gas disks is 

Q = c sK /irGZ, (1) 

where c s is the sound speed, k is the epicyclic frequency at 
which a fluid element oscillates when perturbed from circu- 
lar motion, G is the gravitational constant, and X is the sur- 
face density. In a nearly Keplerian disk, k w the rotational 
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angular speed fl For axisymmetric (ring-like) disturbances, 
disks are stable when Q > 1 (Toomre, 1964). At high Q- 
values, pressure, represented by c s in (1), stabilizes short 
wavelengths, and rotation, represented by k, stabilizes long 
wavelengths. The most unstable wavelength when Q < 1 
is given by A m w 2tt 2 G^/k 2 . 

Modern numerical simulations, beginning with Pa- 
paloizou and Savonije (1991), show that nonaxisymmetric 
disturbances, which grow as multi-armed spirals, become 
unstable for Q < 1.5. Because the instability is both linear 
and dynamic, small perturbations grow exponentially on the 
time scale of a rotation period P rot = 2w/il. The multi-arm 
spiral waves that grow have a predominantly trailing pat- 
tern, and several modes can appear simultaneously (Boss, 
1998a; Laughlin et al, 1998; Nelson et al, 1998; Pickett et 
al., 1998). Although the star does become displaced from 
the system center of mass (Rice et al, 2003a) and one- 
armed structures can occur (see Fig. 1 of Cai et al, 2006), 
one-armed modes do not play the dominant role predicted 
by Adams et al (1989) and Shu et al (1990). 

2.2 Nonlinear Regime 

Numerical simulations (see also Sections 3 and 4) show 
that, as GI's emerge from the linear regime, they may either 
saturate at nonlinear amplitude or fragment the disk. Two 
major effects control or limit the outcome - disk thermo- 
dynamics and nonlinear mode coupling. At this point, the 
disks also develop large surface distortions. 

Disk Thermodynamics. As the spiral waves grow, they 
can steepen into shocks that produce strong localized heat- 
ing (Pickett et al, 1998, 2000a; Nelson et al, 2000). Gas 
is also heated by compression and through net mass trans- 
port due to gravitational torques. The ultimate source of 
GI heating is work done by gravity. What happens next de- 
pends on whether a balance can be reached between heating 
and the loss of disk thermal energy by radiative or convec- 
tive cooling. The notion of a balance of heating and cooling 
in the nonlinear regime was described as early as 1965 by 
Goldreich and Lynden-Bell and has been used as a basis 
for proposing a-treatments for Gl-active disks (Paczynski, 
1978; Lin and Pringle, 1987). For slow to moderate cool- 
ing rates, numerical experiments, such as in Fig. 1, verify 
that thermal self-regulation of GI's can be achieved (Tom- 
ley et al, 1991; Pickett et al, 1998, 2000a, 2003; Nelson 
et al, 2000; Gammie, 2001; Boss, 2003; Rice et al, 2003b; 
Lodato and Rice, 2004, 2005; Mejia et al. 2005; Cai et 
al, 2006). Q then hovers near the instability limit, and the 
nonlinear amplitude is controlled by the cooling rate. 

Nonlinear Mode Coupling. Using second and third- 
order governing equations for spiral modes and comparing 
their results with a full nonlinear hydrodynamics treatment, 
Laughlin et al. (1997, 1998) studied nonlinear mode cou- 
pling in the most detail. Even if only a single mode ini- 
tially emerges from the linear regime, power is quickly dis- 
tributed over modes with a wide variety of wavelengths and 
number of arms, resulting in a self-gravitating turbulence 
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Fig. 1. — Greyscale of effective temperature T e ff in degrees 
Kelvin for a face-on Gl-active disk in an asymptotic state of ther- 
mal self-regulation. This figure is for the Mejia et al. (2005) evo- 
lution of a 0.07 M© disk around a 0.5 Mq star with t coo i — 1 
outer rotation period at 4,500 yr. The frame is 120 AU on a side. 



that permeates the disk. In this gravitoturbulence, gravita- 
tional torques and even Reynold's stresses may be impor- 
tant over a wide range of scales (Nelson et al, 1998; Gam- 
mie, 2001; Lodato and Rice, 2004; Mejia et al, 2005). 

Surface Distortions. As emphasized by Pickett et al. 
(1998, 2000, 2003), the vertical structure of the disk plays a 
crucial role, both for cooling and for essential aspects of the 
dynamics. There appears to be a relationship between GI 
spiral modes and the surface or f-modes of stratified disks 
(Pickett et al, 1996; Lubow and Ogilvie, 1998). As a re- 
sult, except for isothermal disks, GI's tend to have large 
amplitudes at the surface of the disk. Shock heating in the 
GI spirals can also disrupt vertical hydrostatic equilibrium, 
leading to rapid vertical expansions that resemble hydraulic 
jumps (Boley et al, 2005; Boley and Durisen, 2006). The 
resulting spiral corrugations can produce observable effects 
(e.g., masers, Durisen et al, 2001). 

2.3 Heating and Cooling 

Protoplanetary disks are expected to be moderately thin, 
with H/r ~ 0.05 — 0.1, where H is the vertical scale 
height and r is the distance from the star. For hydro- 
static equilibrium in the vertical direction, H w c s /fl. 
The ratio of disk internal energy to disk binding energy 
~ c s 2 /(rO) 2 ~ (H/r) 2 is then < 1%. As growing modes 
become nonlinear, they tap the enormous store of gravita- 
tional energy in the disk. Simulation of the disk energy bud- 
get must be done accurately and include all relevant effects, 
because it is the disk temperature, through c s in equation 
1, that determines whether the disk becomes or remains un- 
stable, once the central mass, which governs most of k, and 
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the disk mass distribution S have been specified. 
2.3.1 Cooling 

There have been three approaches to cooling - make 
simple assumptions about the equation of state (EOS), in- 
clude idealized cooling characterized by a cooling time, or 
treat radiative cooling using realistic opacities. 

EOS. This approach has been used to study mode cou- 
pling (e.g., Laughlin et al. 1998) and to examine disk frag- 
mentation in the limits of isentropic and isothermal behav- 
ior (e.g., Boss, 1998a, 2000; Nelson et al, 1998; Pickett et 
al, 1998, 2003; Mayer et al, 2004). Isothermal evolution 
of a disk, where the disk temperature distribution is held 
fixed in space or when following fluid elements, effectively 
assumes rapid loss of energy produced by shocks and PdV 
work. Isentropic evolution, where specific entropy is held 
fixed instead of temperature, is a more moderate assump- 
tion but is still lossy because it ignores entropy generation 
in shocks (Pickett et al, 1998, 2000a). Due to the energy 
loss, we do not refer to such calculations as adiabatic. Here, 
we restrict adiabatic evolution to mean cases where the fluid 
is treated as an ideal gas with shock heating included via an 
artificial viscosity term in the internal energy equation but 
no radiative cooling. Such calculations are adiabatic in the 
sense that there is no energy loss by the system. Examples 
include a simulation in Pickett et al. (1998) and simulations 
in Mayer et al. (2002, 2004). Mayer et al. use adiabatic 
evolution throughout some simulations, but, in others that 
are started with a locally isothermal EOS, they switch to 
adiabatic evolution as the disk approaches fragmentation. 

Simple Cooling Laws. Better experimental control over 
energy loss is obtained by adopting simple cooling rates per 
unit volume A = e/t coo i, where e is the internal energy per 
unit volume. The t coo i is specified either as a fixed frac- 
tion of the local disk rotation period P ro t, usually by set- 
ting t coo iVl = constant (Gammie, 2001; Rice et al, 2003b; 
Mayer et al, 2004b, 2005) or t coo i = constant everywhere 
(Pickett et al, 2003; Mejia et al, 2005). In the Mayer et 
al. work, the cooling is turned off in dense regions to sim- 
ulate high optical depth. Regardless of t coo i prescription, 
the amplitude of the GI's in the asymptotic state (see Fig. 
1), achieved when heating and cooling are balanced, in- 
creases as t coo i decreases. In addition to elucidating the 
general physics of GI's, such studies address whether GI's 
are intrinsically a local or global phenomenon (Laughlin 
and Rdzyczka 1996; Balbus and Papaloizou, 1999) and 
whether they can be properly modeled by a simple a pre- 
scription. When t coo i is globally constant, the transport in- 
duced by GI's is global with high mass inflow rates (Mejia 
et al, 2005; Michael et al, in preparation); when t coo i^l 
is constant, transport is local, except for thick or very mas- 
sive disks, and the inflow rates are well characterized by a 
constant a (Gammie, 2001; Lodato and Rice, 2004, 2005). 

Radiative Cooling. The published literature on this so far 
comes from only three research groups (Nelson et al, 2000; 
Boss, 2001, 2002b, 2004a; Mejia, 2004; Cai et al, 2006), 
but work by others is in progress. Because Solar System- 
sized disks encompass significant volumes with small and 



large optical depth, this becomes a difficult 3D radiative 
hydrodynamics problem. Techniques will be discussed in 
Section 3.2. For a disk spanning the conventional planet- 
forming region, the opacity is due primarily to dust. Com- 
plications which have to be considered include the dust size 
distribution, its composition, grain growth and settling, and 
the occurrence of fast cooling due to convection. 

In general, the radiative cooling time is dependent on 
temperature T and metallicity Z. Let n r <~ ZT@ r and 
k p ~ ZT^p be the Rosseland and Planck mean opacities, 
respectively, and let r ~ n r H be the vertical optical depth 
to the midplane. For large r, 

tcooi ~ T/Tc f f 4 ~ T- 3 r ~ T~ 2 - 5 +^Z; (2) 

for small r, 

tcooi ~ T/ Kp T 4 ~ T- 3 -^/Z. (3) 

For most temperatures regimes, we expect — 3 < [3 < 2.5, 
so that tcooi increases as T decreases. As Z increases, t coo i 
increases in optically thick regions, but decreases in opti- 
cally thin ones. 
2.3.2 Heating 

In addition to the internal heating caused by GI's through 
shocks, compression, and mass transport, there can be heat- 
ing due to turbulent dissipation (Nelson et al, 2000) and 
other sources of shocks. In addition, a disk may be exposed 
to one or more external radiation fields due to a nearby 
OB star (e.g., Johnstone et al, 1998), an infalling envelope 
(e.g., D'Alessio et al, 1997), or the central star (e.g., Chi- 
ang and Goldreich, 1997). These forms of heat input can be 
comparable to or larger than internal sources of heating and 
can influence Q and the surface boundary conditions. Only 
crude treatments have been done so far for envelope irradi- 
ation (Boss 2001, 2002b; Cai et al, 2006) and for stellar 
irradiation (Mejia, 2004). 

2.4 Fragmentation 

As shown first by Gammie (2001) for local thin-disk cal- 
culations and later confirmed by Rice et al. (2003b) and 
Mejia et al. (2005) in full 3D hydro simulations, disks with 
a fixed tcooi fragment for sufficiently fast cooling, specifi- 
cally when tcooiQ ^ 3, or, equivalently, t coo i < P ro t/2. Fi- 
nite thickness has a slight stabilizing influence (Rice et al, 
2003b; Mayer et al, 2004a). When dealing with realistic ra- 
diative cooling, one cannot apply this simple fragmentation 
criterion to arbitrary initial disk models. One has to apply 
it to the asymptotic phase after nonlinear behavior is well- 
developed (Johnson and Gammie, 2003). Cooling times can 
be much longer in the asymptotic state than they are ini- 
tially (Cai et al, 2006, Mejia et al, in preparation). For 
disks evolved under isothermal conditions, where a simple 
cooling time cannot be defined, local thin-disk calculations 
show fragmentation when Q < 1.4 (Johnson and Gammie, 
2003). This is roughly consistent with results from global 
simulations (e.g., Boss, 2000; Nelson et al, 1998; Pickett et 
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Fig. 2. — Midplane density contours for the isothermal evolution 
of a 0.09 Mq disk around a 1 Mq star. A multi-Jupiter mass 
clump forms near 12 o'clock by 374 years. The frame in the figure 
is 40 AU on a side. The figure is adapted from Boss (2000). 

al., 2000a, 2003; Mayer et al, 2002, 2004a). Fig. 2 shows 
a classic example of a fragmenting disk. 

Although there is agreement on conditions for fragmen- 
tation, two important questions remain. Do real disks ever 
cool fast enough for fragmentation to occur, and do the 
fragments last long enough to contract into permanent pro- 
toplanets before being disrupted by tidal stresses, shear 
stresses, physical collisions, and shocks? 

3. NUMERICAL METHODS 

A full understanding of disk evolution and the planet for- 
mation process cannot easily be obtained using a purely an- 
alytic approach. Although numerical methods are power- 
ful, they have flaws and limitations that must be taken into 
account when interpreting results. Here we describe some 
commonly used numerical techniques and their limitations. 

3.1 Hydrodynamics 

Numerical models have been implemented using one or 
the other of two broad classes of techniques to solve the 
hydrodynamic equations. Each class discretizes the system 
in fundamentally different ways. On one hand, there are 
particle-based simulations using Smoothed Particle Hydro- 
dynamics (SPH) (Benz, 1990; Monaghan, 1992), and, on 
the other, grid-based techniques (e.g., Tohline, 1980; Fryx- 
ell et al, 1991; Stone and Norman, 1992; Boss and Myhill, 
1992; Pickett, 1995). 

SPH uses a collection of particles distributed in space 
to represent the fluid. Each particle is free to move in re- 
sponse to forces acting on it, so that the particle distribu- 
tion changes with the system as it evolves. The particles 
are collisionless, meaning that they do not represent actual 
physical entities, but rather points at which the underlying 
distributions of mass, momentum, and energy are sampled. 
In order to calculate hydrodynamic quantities such as mass 



density or pressure forces, contributions from other parti- 
cles within a specified distance, the smoothing length, are 
weighted according to a smoothing kernel and summed in 
pairwise fashion. Mutual gravitational forces are calculated 
by organizing particles into a tree, where close particles are 
treated more accurately than aggregates on distant branches. 

Grid-based methods use a grid of points, usually fixed in 
space, on which fluid quantities are defined. In the class of 
finite difference schemes, fluxes of mass, momentum, and 
energy between adjacent cells are calculated by taking fi- 
nite differences of the fluid quantities in space. Although 
not commonly used in simulations of GI's, the Piecewise 
Parabolic Method (PPM) of Collela and Woodward (1984) 
represents an example of the class of finite volume schemes. 
For our purposes, an important distinguishing factor is that 
while finite difference and SPH methods may require ar- 
tificial viscosity terms to be added to the equations to en- 
sure numerical stability and produce correct dissipation in 
shocks, PPM does not. 

3.2 Radiative Physics 

In Section 2.3, we describe a number of processes by 
which disks may heat and cool. In this section, we discuss 
various code implementations and their limitations. 

Fixed EOS evolution is computationally efficient be- 
cause it removes the need to solve an equation for the en- 
ergy balance. On the other hand, the gas instantly radiates 
away all heating due to shocks and, for the isothermal case, 
due to compressional heating as well. As a consequence, 
the gas may compress to much higher densities than are re- 
alistic, biasing a simulation towards GI growth and frag- 
mentation even when a physically appropriate temperature 
or entropy scale is used. Although fixed i coo /'s represent 
a clear advance over fixed EOS's, equations 2 and 3 show 
that increasing the temperature, which makes the disk more 
stable, also decreases t coo i. So it is incorrect to view short 
global cooling times as necessarily equivalent to more rapid 
GI growth and fragmentation. In order for fragmentation to 
occur, one needs both a short t coo i and a disk that is cool 
enough to be unstable (e.g., Rafikov, 2005). 

The most physically inclusive simulations to date em- 
ploy radiative transport schemes that allow t coo i to be de- 
termined by disk opacity. Current implementations (Section 
2.3) employ variants of a radiative diffusion approximation 
in regions of medium to high optical depth r, integrated 
from infinity toward the disk midplane. On the other hand, 
radiative losses actually occur from regions where r < 1, 
and so the treatment of the interface between optically thick 
and thin regions strongly influences cooling. Three groups 
have implemented different approaches. 

Nelson et al. (2000) assume that the vertical structure 
of the disk can be defined at each point as an atmosphere 
in thermal equilibrium. In this limit, the interface can 
be defined by the location of the disk photosphere, where 
t = 2/3 (see, e.g., Mihalas, 1977). Cooling at each point 
is then defined as that due to a blackbody with the temper- 
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ature of the photosphere. Boss (2001, 2002b, 2004a, 2005) 
performs a 3D flux-limited radiative diffusion treatment for 
the optically thick disk interior (Bodenheimer et al, 1990), 
coupled to an outer boundary condition where the tempera- 
ture is set to a constant for r < 10, t being measured along 
the radial direction. Meji'a (2004) and Cai et al. (2006) 
use the same radiative diffusion treatment as Boss in their 
disk interior, but they define the interface using r = 2/3, 
measured vertically, above which an optically thin atmo- 
sphere model is self-consistently grafted onto the outward 
flux from the interior. As discussed in Section 4.2, results 
for the three groups differ markedly, indicating that better 
understanding of radiative cooling at the disk surface will 
be required to determine the fate of GI's. 

3.3 Numerical Issues 

The most important limitations facing numerical simula- 
tions are finite computational resources. Simulations have 
a limited duration with a finite number of particles or cells, 
and they must have boundary conditions to describe behav- 
ior outside the region being computed. A simulation must 
distribute grid cells or particles over the interesting parts of 
the system to resolve the relevant physics and avoid errors 
associated with incorrect treatment of the boundaries. Here 
we describe a number of requirements for valid simulations 
and pitfalls to be avoided. 

For growth of GI's, simulations must be able to resolve 
the wavelengths of the instabilities underlying the fragmen- 
tation. Bate and Burke rt (1997) and Truelove et al. (1997) 
each define criteria based on the collapse of a Jeans unsta- 
ble cloud that links a minimum number of grid zones or 
particles to either the physical wavelength or mass associ- 
ated with Jeans collapse. Nelson (2006) notes that a Jeans 
analysis may be less relevant for disk systems because they 
are flattened and rotating rather than homogeneous and in- 
stead proposes a criterion based on the Toomre wavelength 
in disks. Generally, grid-based simulations must resolve the 
appropriate local instability wavelength with a minimum of 
4 to 5 grid zones in each direction, while SPH simulations 
must resolve the local Jeans or Toomre mass with a mini- 
mum of a few hundred particles. 

Resolution of instability wavelengths will be insufficient 
to ensure validity if either the hydrodynamics or gravita- 
tional forces are in error. For example, errors in the hydro- 
dynamics may develop in SPH and finite difference meth- 
ods because a viscous heating term must be added artifi- 
cially to model shock dissipation and, in some cases, to en- 
sure numerical stability. In practice, the magnitude of dissi- 
pation depends in part on cell dimensions rather than just on 
physical properties. Discontinuities may be smeared over 
as many as <~ 10 or more cells, depending on the method. 
Further, Mayer et al. (2004a) have argued that because it 
takes the form of an additional pressure, artificial viscos- 
ity may by itself reduce or eliminate fragmentation. On the 
other hand, artificial viscosity can promote the longevity of 
clumps (see Fig. 3 of Durisen, 2006). 



Gravitational force errors develop in grid simulations 
from at least two sources. First, when Pickett et al. (2003) 
place a small blob within their grid, errors occur in the self- 
gravitation force of the blob that depend on whether the 
cells containing it have the same spacing in each coordi- 
nate dimension. Ideally, grid zones would have comparable 
spacing in all directions, but disks are both thin and radi- 
ally extended. Use of spherical and cylindrical grids tends 
to introduce disparity in grid spacing. Second, Boss (2000) 
shows that maximum densities inside clumps are enhanced 
by orders of magnitude as additional terms in his Poisson 
solver, based on a Y; TO decomposition, are included. SPH 
simulations encounter a different source of error because 
gravitational forces must be softened in order to preserve 
the collisionless nature of the particles. Bate and Burkert 
(1997) and Nelson (2006) each show that large imbalances 
between the gravitational and pressure forces can develop if 
the length scales for each are not identical, possibly induc- 
ing fragmentation in simulations. On the other hand, spa- 
tially and temporally variable softening implies a violation 
of energy conservation. Quantifying errors from sources 
such as insufficiently resolved shock dissipation or gravita- 
tional forces cannot be reliably addressed except by exper- 
imentation. Results of otherwise identical simulations per- 
formed at several resolutions must be compared, and iden- 
tical models must be realized with more than one numeri- 
cal method (as in Section 4.4), so that deficiencies in one 
method can be checked against strengths in another. 

The disks relevant for GI growth extend over several or- 
ders of magnitude in radial range, while GI's may develop 
large amplitudes only over some fraction of that range. 
Computationally affordable simulations therefore require 
both inner and outer radial boundaries, even though the disk 
may spread radially and spiral waves propagate up to or be- 
yond those boundaries. In grid-based simulations, Pickett et 
al. (2000b) demonstrate that numerically induced fragmen- 
tation can occur with incorrect treatment of the boundary. 
Studies of disk evolution must ensure that treatment of the 
boundaries does not produce artificial effects. 

In particle simulations, where there is no requirement 
that a grid be fixed at the beginning of the simulation, 
boundaries are no less a problem. The smoothing in 
SPH requires that the distribution of neighbors over which 
the smoothing occurs be relatively evenly distributed in a 
sphere around each particle for the hydrodynamic quanti- 
ties to be well defined. At currently affordable resolutions 
(<~ 10 5 — 10 6 particles), the smoothing kernel extends over a 
large fraction of a disk scale height, so meeting this require- 
ment is especially challenging. Impact on the outcomes of 
simulations has not yet been quantified. 

4. KEY ISSUES 
4.1 Triggers for GI's 

When disks become unstable, they may either fragment 
or enter a self-regulated phase depending on the cooling 
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Fig. 3. — Face-on density maps for two simulations of interacting M = O.1M protoplanetary disks in binaries with t coo i — 0.5P rot 
viewed face-on. The binary in the left panel has a nearly circular binary orbit with an initial separation of 60 AU and is shown after first 
pericentric passage at 150 yrs (left) and then at 450 yrs (right). Large tidally induced spiral arms are visible at 150 yrs. The right panel 
shows a snapshot at 160 yrs from a simulation starting from an initial orbital separation that is twice as large. In this case, fragmentation 
into permanent clumps occurs after a few disk orbital times. Figures adapted from Mayer et al. (2005). 



time. It is therefore important to know how and when GI's 
may arise in real disks and the physical state of the disk at 
that time. Various mechanisms for triggering GI's are con- 
ceivable, but only a few have yet been studied in any detail. 
Possibilities include - the formation of a massive disk from 
the collapse of a protostellar cloud (e.g., Laughlin and Bo- 
denheimer, 1994; Yorke and Bodenheimer, 1999), clumpy 
infall onto a disk (Boss, 1997, 1998a), cooling of a disk 
from a stable to an unstable state, slow accretion of mass, 
accumulation of mass in a magnetically dead zone, pertur- 
bations by a binary companion, and close encounters with 
other star/disk systems (Boffin et al, 1998; Lin et al, 1998). 
A few of these will be discussed further, with an emphasis 
on some new results on effects of binarity. 

Several authors start their disks with stable or marginally 
stable Q-values and evolve them to instability either by slow 
idealized cooling (e.g., Gammie, 2001; Pickett et al, 2003; 
Mejia et al, 2005) or by more realistic radiative cooling 
(e.g., Johnson and Gammie, 2003; Boss, 2005, 2006; Cai 
et al, 2006). To the extent tested, fragmentation in ideal- 
ized cooling cases are consistent with the Gammie criterion 
(Section 2.4). With radiative cooling, as first pointed out by 
Johnson and Gammie (2003), it is difficult to judge whether 
a disk will fragment when it reaches instability based on its 
initial t coo i. When Mayer et al. (2004a) grow the mass of a 
disk while keeping its temperature constant, dense clumps 
form in a manner similar to clump formation starting from 
an unstable disk. A similar treatment of accretion needs 
to be done using realistic radiative cooling. Simulations 
like these suggest that, in the absence of a strong additional 
source of heating, GI's are unavoidable in protoplanetary 
disks with sufficient mass (<~ 0.1M Q for a ~ 1M Q star). 

A disk evolving primarily due to magnetorotational in- 
stabilites (MRI's) may produce rings of cool gas in the disk 
midplane where the ionization fraction drops sufficiently to 
quell MRI's (Gammie, 1996; Fleming and Stone, 2003). 
Dense rings associated with these magnetically dead zones 



should become gravitationally unstable and may well trig- 
ger a localized onset of GI's. This process might lead to 
disk outbursts related to FU Orionis events (Armitage et al, 
2001) and induce chondrule-forming episodes (Boley et al, 
2005). 

A phase of GI's robust enough to lead to gas giant pro- 
toplanet formation might be achieved through external trig- 
gers, like a binary star companion or a close encounter with 
another protostar and its disk. A few studies have explored 
the effects of binary companions on GI's. Nelson (2000) 
follows the evolution of disks in an equal-mass binary sys- 
tem with a semimajor axis of 50 AU and an eccentricity of 
0.3 and finds that the disks are heated by internal shocks and 
viscous processes to such an extent as to become too hot for 
gas giant planet formation either by disk GI's or by core ac- 
cretion, because volatile ices and organics are vaporized. In 
a comparison of the radiated emission calculated from his 
simulation to those from the L1551 IRS5 system, Nelson 
(2000) finds that the simulation is well below the observed 
system and therefore that the temperatures in the simulation 
are underestimates. He therefore concludes that "planet for- 
mation is unlikely in equal-mass binary systems with a ~ 
50 AU." Currently, over two dozen binary or triple star sys- 
tems have known extrasolar planets, with binary separations 
ranging from <~ 10 AU to <~ 10 3 AU, so some means must 
be found for giant planet formation in binary star systems 
with relatively small semimajor axes. 

Using idealized cooling, Mayer et al. (2005) find that 
the effect of binary companions depends on the mass of 
the disks involved and on the disk cooling rate. For a pair 
of massive disks (M <~ 0.1M Q ), formation of permanent 
clumps can be suppressed as a result of intense heating 
from spiral shocks excited by the tidal perturbation (Fig. 
3 left panel). Clumps do not form in such disks for bi- 
nary orbits having a semimajor axis of ~ 60 AU even when 
tcooi < Prot- The temperatures reached in these disks are 
> 200 K and would vaporize water ice, hampering core ac- 
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cretion, as argued by Nelson (2000). On the other hand, 
pairs of less massive disks (M ~ 0.05M Q ) that would not 
fragment in isolation since they start with Q ~ 2, can pro- 
duce permanent clumps provided that t coo i < P rot . This is 
because the tidal perturbation is weaker in this case (each 
perturber is less massive) and the resulting shock heating is 
thus diminished. Finally, the behavior of such binary sys- 
tems approaches that seen in simulations of isolated disks 
once the semimajor axis grows beyond 100 AU (Fig. 3 right 
panel). 

Calculations by Boss (2006) of the evolution of initially 
marginally gravitationally stable disks show that the pres- 
ence of a binary star companion could help to trigger the 
formation of dense clumps. The most likely explanation 
for the difference in outcomes between the models of Nel- 
son (2000) and Boss (2006) is the relatively short cooling 
times in the latter models (<~ 1 to 2 P ro t, see Boss 2004a) 
compared to the effective cooling time in Nelson (2000) of 
<~ !5P rot at 5 AU, dropping to <~ P rot at 15 AU. Similarly, 
some differences in outcomes between the results of Boss 
(2006) and Mayer et al. (2005) can be expected based on 
different choices of the binary semimajor axes and eccen- 
tricities and differences in the thermodynamics. For exam- 
ple, Mayer et al. (2005) turn off cooling in regions with 
densities higher than 10~ 10 g cm~ 3 to account for high op- 
tical depths. 

Overall, the three different calculations agree that exci- 
tation or suppression of fragmentation by a binary compan- 
ion depends sensitively on the balance between compres- 
sional/shock heating and cooling. This balance appears to 
depend on the mass of the disks involved. Interestingly, 
lighter disks are more likely to fragment in binary systems 
according to both Mayer et al. (2005) and Boss (2006). 

4.2 Disk Thermodynamics 

As discussed in Sections 2.2 and 4.1, heating and cool- 
ing are perhaps the most important processes affecting the 
growth and fate of GI's. Thermal regulation in the nonlinear 
regime leads naturally to systems near their stability limit 
where temporary imbalances in one heating or cooling term 
lead to a proportionate increase in a balancing term. For 
fragmentation to occur, a disk must cool quickly enough, 
or fail to be heated for long enough, to upset this self- 
regulation. A complete model of the energy balance that 
includes all relevant processes in a time-dependent manner 
is beyond the capabilities of the current generation of mod- 
els. It requires knowledge of all the following - external ra- 
diation sources and their influence on the disk at each loca- 
tion, the energy lose rate of the disk due to radiative cooling, 
dynamical processes that generate thermal energy through 
viscosity or shocks, and a detailed equation of state to deter- 
mine how much heating any of those dynamical processes 
generate. Recent progress towards understanding disk evo- 
lution has focused on the more limited goals of quantifying 
the sensitivity of results to various processes in isolation. 

In a thin, steady state a-disk, the heating and cooling 



times are the same and take a value (Pringle, 1981; Gam- 
mie, 2001): 

t coo i= I [ 7 (7-l)afi] _1 - (4) 

For a ~ 10~ 2 and 7 = 1.4, equation 4 gives <~ 12P rot . 
This is a crude upper limit on the actual time scale required 
to change the disk thermodynamic state. External radia- 
tive heating from the star and any remaining circumstellar 
material can contribute a large fraction of the total heating 
(D'Alessio et al, 1998; Nelson et al, 2000), as will any 
internal heating due to globally generated dynamical insta- 
bilities that produce shocks. Each of these processes actu- 
ally makes the disk more stable by heating it, but, as a con- 
sequence, dynamical evolution slows until the disk gains 
enough mass to become unstable again. The marginally sta- 
ble state will then be precariously held because the higher 
temperatures mean that all of the heating and cooling time 
scales, i.e., the times required to remove or replace all the 
disk thermal energy, are short (equations 2 and 3). When 
the times are short, any disruption of the contribution from 
a single source may be able to change the thermodynamic 
state drastically within only a few orbits, perhaps beyond 
the point where balance can be restored. 

A number of models (Section 2.3) have used fixed EOS 
evolution instead of a full solution of an energy equation 
to explore disk evolution. A fixed EOS is equivalant to 
specifying the outcomes of all heating and cooling events 
that may occur during the evolution, short-circuiting ther- 
mal feedback. If, for example, the temperature or entropy 
is set much too high or too low, a simulation may predict 
either that no GI's develop in a system, or that they in- 
evitably develop and produce fragmentation, respectively. 
Despite this limitation, fixed EOS's have been useful to de- 
lineate approximate boundaries for regions of marginal sta- 
bility. Since the thermal state is fixed, disk stability (as 
quantified by equation 1) is essentially determined by the 
disk's mass and spatial dimensions, though its surface den- 
sity. Marginal stability occurs generally at Q w 1.2 to 1.5 
for locally isentropic evolutions, with a tendency for higher 
<5's being required to ensure stability with softer EOS's 
(i.e., with lower 7 values) (Boss 1998a; Nelson et al, 1998; 
Pickett et al, 1998, 2000a; Mayer et al, 2004a). At tem- 
peratures appropriate for observed systems (e.g., Beckwith 
etal, 1990), these Q values correspond to disks more mas- 
sive than ~ 0.1M* or surface densities E > 10 3 gm/cm 2 . 

As with their fixed EOS cousins, models with fixed t coo i 
can quantify boundaries at which fragmentation may set in. 
They represent a clear advance over fixed EOS evolution 
by allowing thermal energy generated by shocks or com- 
pression to be retained temporarily, and thereby enabling 
the disk's natural thermal regulation mechanisms to deter- 
mine the evolution. Models that employ fixed cooling times 
can address the question of how violently the disk's ther- 
mal regulation mechanisms must be disrupted before they 
can no longer return the system to balance. An example 
of the value of fixed t coo i calculations is the fragmentation 
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criterion t coo i < (see Section 2.4). 

The angular momentum transport associated with disk 
self-gravity is a consequence of the gravitational torques 
induced by GI spirals (e.g., Larson, 1984). The viscous 
a parameter is actually a measure of that stress normal- 
ized by the local disk pressure. As shown in equation 4 
and reversing the positions of t coo i and a, the stress in a 
self-gravitating disk depends on the cooling time and on 
the equation of state through the specific heat ratio. As 
long as the dimensionless scale height is H < 0.1, global 
simulations by Lodato and Rice (2004) with t coo iVl = con- 
stant confirm Gammie's assumption that transport due to 
disk self-gravity can be modeled as a local phenomenon 
and that equation 4 is accurate. Gammie (2001) and Rice et 
al. (2005) show that there is a maximum stress that can be 
supplied by such a quasi-steady, self-gravitating disk. Frag- 
mentation occurs if the stress required to keep the disk in 
a quasi-steady state exceeds this maximum value. The re- 
lationship between the stress and the specific heat ratio, 7, 
results in the cooling time required for fragmentation in- 
creasing as 7 decreases. For 7 = 7/5, the cooling time 
below which fragmentation occurs may be more like 2P rot , 
not the 3/51 ~ P r ot/2 obtained for 7 = 2 (Gammie, 2001; 
Mayer et al., 2004b; Rice et al, 2005). 

Important sources of stress and heating in the disk, 
that lie outside the framework of Gammie's local analysis, 
are global gravitational torques due to low-order GI spiral 
modes. There are two ways this can happen - a geomet- 
rically thick massive disk {Lodato and Rice, 2005) and a 
fixed global t coo i = constant (Mejia et al., 2005). Disks 
then initially produce large-amplitude spirals, resulting in a 
transient burst of global mass and angular momentum re- 
distribution. For t coo i = constant and moderate masses, the 
disks then settle down to a self -regulated asymptotic state 
but with gravitational stresses significantly higher than pre- 
dicted by equation 4 (Michael et al. , in preparation). For the 
very massive t coo i^ = constant disks, recurrent episodic re- 
distributions occur. In all these cases, the heating in spiral 
shocks is spatially and temporally very inhomogeneous, as 
are fluctuations in all thermodynamic variables and the ve- 
locity field. 

The most accurate method to determine the internal ther- 
modynamics of the disk is to couple the equations of radia- 
tive transport to the hydrodynamics directly. All heating or 
cooling due to radiation will then be properly defined by 
the disk opacity, which depends on local conditions. This 
is important because some fraction of the internal heating 
will be highly inhomogeneous, occurring predominantly in 
compressions and shocks as gas enters a high density spiral 
structure, or at high altitudes where waves from the interior 
are refracted and steepen into shocks (Pickett et al. 2000a) 
and where disks may be irradiated (Mejia, 2004; Cai et al., 
2006). Temperatures and the t coo is that depend on them 
will then be neither simple functions of radius, nor a sin- 
gle globally defined value. Depending on whether the lo- 
cal cooling time of the gas inside the high density spiral 
structure is short enough, fragmentation will be more or less 



likely, and additional hydrodynamic processes such as con- 
vection may become active if large enough gradients can be 
generated. 

Indeed, recent simulations of Boss (2002a, 2004a) sug- 
gest that vertical convection is active in disks when radia- 
tive transfer is included, as expected for high r according to 
Ruden and Pollack (1991). This is important because con- 
vection will keep the upper layers of the disk hot, at the 
expense of the dense interior, so that radiative cooling is 
more efficient and fragmentation is enhanced. The results 
have not yet been confirmed by other work and therefore re- 
main somewhat controversial. Simulations by Mejia (2004) 
and Cai et al. (2006) are most similar to those of Boss 
and could have developed convection sufficient to induce 
fragmentation, but none seems to occur. No fragmentation 
occurs in Nelson et al. (2000) either, where convection is 
implicitly assumed to be efficient through their assumption 
that the entropy of each vertical column is constant. Re- 
cent re-analysis of their results reveals t coo i ~ 3 to 10 P ro t, 
depending on radius, which is too long to allow fragmen- 
tation. These i coo z's are in agreement with those seen by 
Cai et al. (2006) and by Mejia et al. (in preparation) for 
solar metallicity. The Nelson et al. results are also interest- 
ing because their comparison of the radiated output to SEDs 
observed for real systems demonstrates that substantial ad- 
ditional heating beyond that supplied by GI's is required to 
reproduce the observations, perhaps further inhibiting frag- 
mentation in their models. However, using the same tem- 
perature distribution between 1 and 10 AU now used in 
Boss's GI models, combined with temperatures outside this 
region taken from models by Adams et al. (1988), Boss and 
Yorke (1996) are able to reproduce the SED of the T Tauri 
system. It is unclear at present why their results differ from 
those of Nelson et al. (2000). 

The origins of the differences between the three studies 
are uncertain, but possibilities include differences of both 
numerical and physical origin. The boundary treatment at 
the optically thick/thin interface is different in each case 
(see Section 3.2), influencing the efficiency of cooling, as 
are the numerical methods and resolutions. Boss and the 
Cai/Mejia group each use 3D grid codes but with spherical 
and cylindrical grids respectively, and each with a different 
distribution of grid zones, while Nelson et al. use a 2D SPH 
code. Perhaps significantly, Cai/Mejia assume their ideal 
gas has 7 = 5/3 while Boss adopts an EOS that includes ro- 
tational and vibrational states of hydrogen, so that 7 m 7/5 
for typical disk conditions. It is possible that differences in 
the current results may be explained if the same sensitivi- 
ties to 7 seen in fixed EOS and fixed cooling simulations 
also hold when radiative transfer is included. Boss and Cai 
(in preparation) are now conducting direct comparison cal- 
culations to isolate the cause of their differences. The pre- 
liminary indication is that the radiative boundary conditions 
may be the critical factor. 

Discrepant results for radiatively cooled models should 
not overshadow the qualitative agreement reached about the 
relationship between disk thermodynamics and fragmenta- 
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tion. If the marginally unstable state of a self-regulated 
disk is upset quickly enough by an increase in cooling or 
decrease in heating, the disk may fragment. What is still 
very unclear is whether such conditions can develop in real 
planet-forming disks. It is key to develop a full 3D por- 
trait of the disk surface, so that radiative heating and cool- 
ing sources may be included self-consistently in numerical 
models. Important heating sources will include the enve- 
lope, the central star, neighboring stars, and self-heating 
from other parts of the disk, all of which will be sensitive 
to shadowing caused by corrugations in the disk surface 
that develop and change with time due to the GI's them- 
selves. Preliminary studies of 3D disk structure (Boley and 
Durisen, 2006) demonstrate that vertical distortions, analo- 
gous to hydraulic jumps, will in fact develop (see also Pick- 
ett et al, 2003). If these corrugations are sufficient to cause 
portions of the disk to be shadowed, locally rapid cooling 
may occur in the shadowed region, perhaps inducing frag- 
mentation. 

An implicit assumption of the discussion above is that 
the opacity is well known. In fact, it is not. The dominant 
source of opacity is dust, whose size distribution, compo- 
sition, and spatial distribution will vary with time (Cuzzi et 
al, 2001; Klahr, 2003, see also Section 5 below), causing 
the opacity to vary as a result. So far, no models of GI evo- 
lution have included effects from any of these processes, ex- 
cept that Nelson et al. model dust destruction while Cai and 
Mejfa consider opacity due to large grains. Possible con- 
sequences are a misidentification of the disk photospheric 
surface if dust grains settle towards the midplane, or incor- 
rect radiative transfer rates in optically thick regions if the 
opacities themselves are in error. 

4.3 Orbital Survival of Clumps 

Once dense clumps form in a gravitationally unstable 
disk, the question becomes one of survival: Are they tran- 
sient structures or permanent precursors of giant planets? 
Long-term evolution of simulations that develop clumps is 
difficult because it requires careful consideration of not only 
the large-scale dynamical processes that dominate forma- 
tion but also physical processes that exert small influences 
over long time scales (e.g., migration and transport due to 
viscosity). It also requires that boundary conditions be han- 
dled gracefully in cases where a clump or the disk itself tries 
to move outside the original computational volume. 

On a more practical level, the extreme computational 
cost of performing such calculations limits the time over 
which systems may be simulated. As a dense clump forms, 
the temperatures, densities, and fluid velocities within it all 
increase. As a result, time steps, limited by the Courant con- 
dition, can decrease to as little as minutes or hours as the 
simulation attempts to resolve the clump's internal struc- 
ture. So far only relatively short integration times of up to 
a fewx 10 3 yrs have been possible. Here, we will focus on 
the results of simulations and refer the reader to the chap- 
ters by Papaloizou et al. and Levison et al. for discussions 



of longer-term interactions. 

In the simplest picture of protoplanet formation via GI's, 
structures are assumed to evolve along a continuum of 
states that are progressively more susceptible to fragmen- 
tation, presumably ending in one or more bound objects 
which eventually become pro toplanets. Pickett et al. (1998, 
2000a, 2003) and Mejfa et al. (2005) simulate initially 
smooth disks subject to growth of instabilities and, in- 
deed, find growth of large-amplitude spiral structures that 
later fragment into arclets or clumps. Instead of growing 
more and more bound, however, these dense structures are 
sheared apart by the background flow within an orbit or less, 
especially when shock heating is included via an artificial 
viscosity. This suggests that a detailed understanding of the 
thermodynamics inside and outside the fragments is critical 
for understanding whether fragmentation results in perma- 
nently bound objects. 

Assuming that permanently bound objects do form, two 
additional questions emerge. First, how do they accrete 
mass and how much do they accrete? Second, how are they 
influenced by the remaining disk material? Recently, Mayer 
et al. (2002, 2004a) and Lufkin et al. (2004) have used 
SPH calculations to follow the formation and evolution of 
clumps in simulations covering up to 50 orbits (roughly 
600 yrs), and Mayer et al. (in preparation) are extending 
these calculations to several thousand years. They find that, 
when a locally isothermal EOS is used well past initial frag- 
mentation, clumps grow to ~ 10Mj within a few hundred 
years. On the other hand, in simulations using an ideal gas 
EOS plus bulk viscosity, accretion rates are much lower 
(< 10~ 6 M Q /yr), and clumps do not grow to more than 
a few Mj or <~ 1% of the disk mass. The assumed ther- 
modynamic treatment has important effects not only on the 
survival of clumps, but also on their growth. 

Nelson and Benz (2003), using a grid-based code and 
starting from a 0.3Mj seed planet, show that accretion rates 
this fast are unphysically high because the newly accreted 
gas cannot cool fast enough, even with the help of convec- 
tion, unless some localized dynamical instability is present 
in the clump's envelope. So, the growth rate of an initially 
small protoplanet may be limited by its ability to accept 
additional matter rather than the disk's ability to supply it. 
They note (see also Lin and Papaloizou, 1993; Bryden et al, 
1999; Kley, 1999; Lubow et al, 1999; Nelson et al, 2000) 
that the accretion process after formation is self-limiting at 
a mass comparable to the largest planet masses yet discov- 
ered (see the chapter by Udry et al). 

Fig. 4 shows one of the extended Mayer et al. sim- 
ulations, containing two clumps in one disk realized with 
2 x 10 5 particles, and run for about 5,000 years (almost 200 
orbits at 10 AU). There is little hint of inward orbital mi- 
gration over a few thousand year time scale. Instead, both 
clumps appear to migrate slowly outward. Boss (2005) uses 
sink particles ("virtual planets") to follow a clumpy disk for 
about 1,000 years. He also finds that the clumps do not mi- 
grate rapidly. In both works, the total simulation times are 
quite short compared to the disk lifetime and so are only 
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Fig. 4. — The orbital evolution of two clumps (right) formed in a massive, growing protoplanetary disk simulation described in Mayer 
et al. (2004). A face-on view of the system after 2264 years of evolution is shown on the left, using a color coded density map (the box 
is 38 AU on a side). In the right panel, the orbital evolution of the two clumps is shown. Overal, both clumps migrate outward. 



suggestive of the longer-term fate of the objects. Neverthe- 
less, the results are important, because they illustrate short- 
comings in current analytic models of migration. 

Although migration theory is now extremely well devel- 
oped (see the chapter by Papaloizou et al), predictions for 
migration at the earliest phases of protoplanet formation by 
GI's are difficult to make, because many of the assumptions 
on which the theory is based are not well satisfied. More 
than one protoplanet may form in the same disk, they may 
form with masses larger than linear theory can accommo- 
date, and they may be significantly extended rather than the 
point masses assumed by theory. If the disk remains mas- 
sive, it may also undergo gravitoturbulence that changes the 
disk's mass distribution on a short enough time scale to call 
into question the resonance approximations in the theory. If 
applicable in the context of these limitations, recent investi- 
gations into the character of corotation resonances (see the 
chapter by Papaloizou et al.) and vortex excitation (Roller 
et al., 2003) in the corotation region may be of particular 
interest, because a natural consequence of these processes 
is significant mass transport across the clump's orbit and re- 
duced inward migration, which is in fact seen in the above 
simulations. 

4.4 Comparison Test Cases 

Disk instability has been studied so far with various 
types of grid codes and SPH codes that have different rel- 
ative strengths and weaknesses (Section 3). Whether dif- 
ferent numerical techniques find comparable results with 
nearly identical assumptions is not yet known, although 
some comparative studies have been attempted (Nelson et 
al., 1998). Several aspects of GI behavior can be highly 
dependent on code type. For example, SPH codes require 
artificial viscosity to handle shocks such as those occurring 
along spiral arms. Numerical viscosity can smooth out the 
velocity field in overdense regions, possibly inhibiting col- 
lapse (Mayer et al, 2004a) but, at the same time, possibly 



increasing clump longevity if clumps form (see Fig. 3 of 
Durisen, 2006). Gravity solvers that are both accurate and 
fast are a robust feature of SPH codes, while gravity solvers 
in grid codes can under-resolve the local self-gravity of the 
gas (Pickett et al, 2003). Both types of codes can lead to 
spurious fragmentation or suppress it when a force imbal- 
ance between pressure and gravity results at scales com- 
parable to the local Jeans or Toomre length due to lack of 
resolution (Truelove et al, 1997; Bate and Burkert, 1997; 
Nelson, 2006). 

Another major code difference is in the set up of ini- 
tial conditions. Although both Eulerian grid-based and La- 
grangian particle-based techniques represent an approxima- 
tion to the continuum fluid limit, noise levels due to dis- 
creteness are typically higher in SPH simulations. Initial 
perturbations are often applied in grid-based simulations to 
seed GI's (either random or specific modes or both, e.g., 
Boss, 1998a), but are not required in SPH simulations, be- 
cause they already have built-in Poissonian noise at the level 
of v^V /N or more, where N is the number of particles. In 
addition, the SPH calculation of hydrodynamic variables in- 
troduces small scale noise at the level of l/N neigh , where 
Nneigh is the number of neighboring particles contained 
in one smoothing kernel. Grid-based simulations require 
boundary conditions which restrict the dynamic range of 
the simulations. For example, clumps may reach the edge 
of a computational volume after only a limited number of 
orbits (Boss, 1998a, 2000; Pickett et al, 2000a). Cartesian 
grids can lead to artificial diffusion of angular momentum 
in a disk, a problem that can be avoided using a cylindrical 
grid (Pickett et al, 2000a) or spherical grid (Boss and My - 
hill, 1992). Myhill & Boss (1993) find good agreement be- 
tween spherical and Cartesian grid results for a nonisother- 
mal rotating protostellar collapse problem, but evolution of 
a nearly equilibrium disk over many orbits in a Cartesian 
grid is probably still a challenge. 

In order to understand how well different numerical tech- 
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Fig. 5. — Equatorial slice density maps of the disk in the test runs after about 100 yrs of evolution. The initial disk is 20 AU in diameter. 
From top left to bottom right are the results from GASOLINE and GADGET2 (both SPH codes), from the Indiana cylindrical-grid code, 
and from the AMR Cartesian-grid code FLASH. The SPH codes adopt the shear-reduced artificial viscosity of Balsam (1995). 



niques can agree on the outcome of GFs, different codes 
need to run the same initial conditions. This is being done 
in a large, on-going code-comparison project that involves 
eight different codes, both grid-based and SPH. Among 
the grid codes, there are several adaptive mesh refinement 
(AMR) schemes. The comparison is part of a larger ef- 
fort involving several areas of computational astrophysics 
(http://krone.physik.unizh.ch/^moore/wengen/tests.html). 
The system chosen for the comparison is a uniform temper- 
ature, massive, and initially very unstable disk with a diam- 
eter of about 20 AU. The disk is evolved isothermally and 
has a Q profile that decreases outward, reaching a minimum 
value ~ 1 at the disk edge. The disk model is created using 
a particle representation by letting its mass grow slowly, as 
described in Mayer et al. (2004a). This distribution is then 
interpolated onto the various grids. 

Here we present the preliminary results of the code com- 
parisons from four codes - two SPH codes called GASO- 
LINE {Wadsley et al, 2004) and GADGET2 (Springel et 
al, 2001; Springel, 2005), the Indiana University code with 
a fixed cylindrical grid (Pickett, 1995; Mejia, 2004), and the 
Cartesian AMR code caUed FLASH (Fryxell et al, 2000). 
Readers should consult the published literature for detailed 
descriptions, but we briefly enumerate some basic features. 
FLASH uses a PPM-based Riemann solver on a Cartesian 
grid with directional splitting to solve the Euler equations, 
and it uses an iterative multi-grid Poisson solver for gravity. 



Both GASOLINE and GADGET2 solve the Euler equations 
using SPH and solve gravity using a treecode, a binary tree 
in the case of GASOLINE and an oct-tree in the case of 
GADGET2. Gravitational forces from individual particles 
are smoothed using a spline kernel softening, and they both 
adopt the Balsara (1995) artificial viscosity that minimizes 
shear forces on large scales. The Indiana code is a finite dif- 
ference grid-based code which solves the equations of hy- 
drodynamics using the Van Leer method. Poisson's equa- 
tion is solved at the end of each hydrodynamic step by a 
Fourier transform of the density in the azimuthal direction, 
direct solution by cyclic reduction of the transform in (r,z), 
and a transform back to real space (Tohline, 1980). The 
code's Von Neumann-Richtmeyer artificial bulk viscosity is 
not used for isothermal evolutions. 

The two SPH codes are run with fixed gravitational soft- 
ening, and the local Jeans length (see Bate and Burkert, 
1997) before and after clump formation is well resolved. 
Runs with adaptive gravitational softening will soon be in- 
cluded in the comparison. Here we show the results of the 
runs whose initial conditions were generated from the 8 x 
10 5 particles setup, which was mapped onto a 512x512x52 
Cartesian grid for FLASH and onto a 512x1024x64 (r,cf>,z) 
cylindrical grid for the Indiana code. Comparable resolu- 
tion (cells for grids or gravity softening for SPH runs) is 
available initially in the outer part of the disk, where the 
Q parameter reaches its minimum. In the GASOLINE and 
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GADGET2 runs, the maximum spatial resolution is set by 
the gravitational softening at 0.12 AU. Below this scale, 
gravity is essentially suppressed. The FLASH run has a 
initial resolution of 0.12 AU at 10 AU, comparable with 
the SPH runs. The Indiana code has the same resolution as 
FLASH in the radial direction but has a higher azimuthal 
resolution of 0.06 AU at 10 AU. 

As it can be seen from Fig. 5, the level of agreement 
between the runs is satisfactory, although significant differ- 
ences are noticeable. More clumps are seen in the Indi- 
ana code simulation. On the other end, clumps have simi- 
lar densities in FLASH and GASOLINE, while they appear 
more fluffy in the Indiana code than in the other three. The 
causes are probably different gravity solvers and the non- 
adaptive nature of the Indiana code. Even within a single 
category of code, SPH or grid-based, different types of vis- 
cosity, both artificial and numerical, might be more or less 
diffusive and affect the formation and survival of clumps. In 
fact, tests show that more fragments are present in SPH runs 
with shear-reduced artificial viscosity than with full shear 
viscosity. 

Although still in an early stage, the code comparison has 
already produced one important result, namely that, once 
favorable conditions exist, widespread fragmentation is ob- 
tained in high-resolution simulations using any of the stan- 
dard numerical techniques. On the other hand, the differ- 
ences already noticed require further understanding and will 
be addressed in a forthcoming paper (Mayer et al, in prepa- 
ration). Although researchers now agree on conditions for 
disk fragmentation, no consensus yet exists about whether 
or where real disks fragment or how long fragments really 
persist. Answers to these questions require advances over 
current techniques for treating radiative physics and com- 
pact structures in global simulations. 

5. INTERACTIONS WITH SOLIDS 

The standard model for the formation of giant gaseous 
planets involves the initial growth of a rocky core that, when 
sufficiently massive, accretes a gaseous envelope (Boden- 
heimer and Pollack, 1986; Pollack et al, 1996). In this 
scenario, the solid particles in the disk must first grow from 
micron-sized dust grains to kilometer-sized planetesimals 
that then coagulate to form the rocky core. 

In a standard protoplanetary disk, the gas pressure near 
the disk midplane will generally decrease with increasing 
radius resulting in an outward pressure gradient that causes 
the gas to orbit with sub-Keplerian velocities. The solid 
particles, on the other hand, do not feel the gas pressure and 
orbit with Keplerian velocities. This velocity difference re- 
sults in a drag force that generally causes the solid particles 
to lose angular momentum and to spiral inward toward the 
central star with a radial drift velocity that depends on the 
particle size (Weidenschilling, 1977). 

While this differential radial drift can mix together par- 
ticles of different size and allow large grains to grow by 
sweeping up smaller grains (Weidenschilling and Cuzzi, 



1993), it also introduces a potential problem. Depending 
on the actual disk properties, the inward radial velocity for 
particles with sizes between 1 cm and 1 m can be as high 
as 10 4 cm s _1 (Weidenschilling, 1977), so that these parti- 
cles could easily migrate into the central star before becom- 
ing large enough to decouple from the disk gas. If these 
particles do indeed have short residence times in the disk, 
it is difficult to envisage how they can grow to form the 
larger kilometer-sized planetesimals which are required for 
the subsequent formation of the planetary cores. 

The above situation is only strictly valid in smooth, lam- 
inar disks with gas pressures that decrease monotonically 
with increasing radius. If there are any regions in the disk 
that have local pressure enhancements, the situation can be 
very different. In the vicinity of a pressure enhancement, 
the gas velocity can be either super- and sub-Keplerian de- 
pending on the local gas pressure gradient. The drag force 
can then cause solid particles to drift outwards or inwards, 
respectively (Haghighipour and Boss, 2003a,b). The net 
effect is that the solid particles should drift towards pres- 
sure maxima. A related idea is that a baroclinic instability 
could lead to the production of long-lived, coherent vortices 
(Klahr and Bodenheimer, 2003) and that solid particles 
would drift towards the center of the vortex where the en- 
hanced concentration could lead to accelerated grain growth 
(Klahr and Henning, 1997). The existence of such vortices 
is, however, uncertain (Johnson and Gammie, 2006). 




Fig. 6. — Surface density structure of particles embedded in a 
self-gravitating gas disk, a) The left-hand panel shows that the 
distribution of 10 m radius particles is similar to that of the gas 
disk, because these particles are not influenced strongly by gas 
drag, b) The right-hand panel illustrates that 50 cm particles are 
strongly influenced by gas drag and become concentrated into the 
Gl-spirals with density enhancements of an order of magnitude or 
more. Figures adapted from Rice et al. (2004). 

An analogous process could occur in a self-gravitating 
disk, where structures formed by GI activity, such as the 
centers of the spiral arms, are pressure and density maxima. 
In such a case, drag force results in solid particles drifting 
towards the centers of these structures, with the most sig- 
nificant effect occurring for those particles that would, in a 
smooth, laminar disk, have the largest inward radial veloci- 
ties. 

If disks around very young protostars do indeed undergo 
a self-gravitating phase, then we would expect the resulting 
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spiral structures to influence the evolution of the solid par- 
ticles in the disk (Haghighipour and Boss, 2003a,b). A GI- 
active disk will also transport dust grains small enough to 
remain tied to the gas across distances of many AU's in only 
1,000 yrs or so (Boss, 2004b), a potentially important pro- 
cess for explaining the components of primitive meteorites 
(see the chapter by Alexander et al). Boley and Durisen 
(2006) show that, in only one pass of a spiral shock, hy- 
draulic jumps induced by shock heating can mix gas and en- 
trained dust radially and vertically over length-scales ~ H 
through the generation of huge breaking waves. The pres- 
ence of chondrules in primitive chondritic meteorites is cir- 
cumstantial evidence that the Solar Nebula experienced a 
self-gravitating phase in which spiral shock waves provided 
the flash heating required to explain their existence (Boss 
and Durisen, 2005a,b; Boley et al, 2005). 

To test how a self-gravitating phase in a protostellar disk 
influences the evolution of embedded particles, Rice et al. 
(2004) perform 3D self-gravitating disk simulations that in- 
clude particles evolved under the influence of both disk self- 
gravity and gas drag. In their simulations, they consider 
both 10 m particles, which, for the chosen disk parameters, 
are only weakly coupled to the gas, and 50 cm particles that 
are significantly influenced by the gas drag. Fig. 6a shows 
the surface density structure of the 10 m particles one outer 
rotation period after they were introduced into the gas disk. 
The structure in the particle disk matches closely that of the 
gas disk (not shown) showing that these particles are influ- 
enced by the gravitational force of the gas disk, but not so 
strongly influenced by gas drag. Fig. 6b shows the sur- 
face density structure of 50 cm particles at the same epoch. 
Particles of this size are influenced by gas drag and Fig. 
6b shows that, compared to the 10 m particles, these parti- 
cles become strongly concentrated into the Gl-induced spi- 
ral structures. 

The ability of solid particles to become concentrated in 
the center of Gl-induced structures suggests that, even if gi- 
ant planets do not form directly via GI's, a self-gravitating 
phase may still play an important role in giant planet forma- 
tion. The solid particles may achieve densities that could 
accelerate grain growth either through an enhanced colli- 
sion rate or through direct gravitational collapse of the par- 
ticle sub-disk (Youdin and Shu, 2002). Durisen et al. (2005) 
also note that dense rings can be formed near the boundaries 
between Gl-active and inactive regions of a disk (e.g., the 
central disk in Fig. 1). Such rings are ideal sites for the con- 
centration of solid particles by gas drag, possibly leading 
to accelerated growth of planetary embryos. Even if pro- 
cesses like these do not contribute directly to planetesimal 
growth, GI's may act to prevent the loss of solids by migra- 
tion toward the proto-Sun. The complex and time-variable 
structure of GI activity should increase the residence time 
of solids in the disk and potentially give them enough time 
to become sufficiently massive to decouple from the disk 
gas. 



6. PLANET FORMATION 

The relatively high frequency (~ 10%) of solar-type stars 
with giant planets that have orbital periods less than a few 
years suggests that longer-period planets may be quite fre- 
quent. Perhaps ~ 12 to 25% of G dwarfs may have gas 
giants orbiting within ~10 AU. If so, gas giant planet for- 
mation must be a fairly efficient process. Because roughly 
half of protoplanetary disks disappear within 3 Myr or less 
(Bally et al, 1998; Haisch et al, 2001 ; Briceno et al, 2001 ; 
Eisner and Carpenter, 2003), core accretion may not be 
able to produce a high frequency of gas giants. There is also 
now strong theoretical (Yorke and Bodenheimer, 1999) and 
observational (Osorio et al, 2003; Rodriguez et al, 2005; 
Eisner et al, 2005) evidence that disks around very young 
protostars should indeed be sufficiently massive to experi- 
ence GI's. Rodriguez et al. (2005) show a 7 mm VLA 
image of a disk around a Class protostar that may have a 
mass half that of the central star. 

Hybrid scenarios may help remove the bottleneck by 
concentrating meter-sized solids, but it is not clear that they 
can shorten the overall time scale for core accretion, which 
is limited by the time needed for the growth of 10 M ffi cores 
and for accretion of a large gaseous envelope. Durisen et 
al. (2005) suggest that the latter might be possible in dense 
rings, but detailed calculations of core growth or envelope 
accretion in the environment of a dense ring do not now 
exist. Disk instability, on the other hand, has no problem 
forming gas giants rapidly in even the shortest-lived pro- 
toplanetary disk. Most stars form in regions of high-mass 
star formation (Lada and Lada, 2003) where disk lifetimes 
should be the shortest due to loss of outer disk gas by UV 
irradiation. 

There is currently disagreement about whether GI's are 
stronger in low-metallicity systems (Cai et al, 2006) or 
whether their strength is relatively insensitive to the opacity 
of the disk (Boss, 2002a). In either case, if disk instability 
is correct, we would expect that even low-metallicity stars 
could host gas giant planets. The growth of cores in the 
core accretion mechanism is hastened by higher metallicity 
through the increase in surface density of solids (Pollack et 
al, 1996), although the increased envelope opacity, which 
slows the collapse of the atmosphere, works in the other di- 
rection (Podolak, 2003). The recent observation of a Saturn 
mass object, orbiting the metal-rich star HD 149026, with 
a core mass equal to approximately half the planet's mass 
(Sato et al, 2005) has been suggested as a strong confir- 
mation of the core accretion model. It has, however, yet to 
be shown that the core accretion model can produce a core 
with such a relatively large mass. If this core was produced 
by core accretion, it seems that it never achieved a runaway 
growth of its envelope; yet, in the case of Jupiter, the core 
accretion scenario requires efficient accumulation of a mas- 
sive envelope around a relatively low-mass core. 

The correlation of short-period gas giants with high 
metallicity stars is often interpreted as strong evidence in 
favor of core accretion (Laws et al, 2003; Fischer et al, 
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2004; Santos et al, 2004). The Santos et al. (2004) anal- 
ysis, however, shows that even the stars with the lowest 
metallicities have detectable planets with a frequency com- 
parable to or higher than that of the stars with intermedi- 
ate metallicities. Rice et al. (2003c) have shown that the 
metallicity distribution of systems with at least one massive 
planet (M p i > 5Mj up ) on an eccentric orbits of moderate 
semi-major axis does not have the same metal-rich nature 
as the full sample of extrasolar planetary systems. Some of 
the metallicity correlation can be explained by the observa- 
tional bias of the spectroscopic method in favor of detecting 
planets orbiting stars with strong metallic absorption lines. 
The residual velocity jitter typically increases from a few 
m/s for solar metallicity to 5 - 16 m/s for stars with 1/4 
the solar metallicity or less. In terms of extrasolar planet 
search space, this could account for as much as a factor 
of two difference in the total number of planets detected 
by spectroscopy. A spectroscopic search of 98 stars in the 
Hyades cluster, with a metallicity 35% greater than solar, 
found nothing, whereas about 10 hot Jupiters should have 
been found, assuming the same frequency as in the solar 
neighborhood (Paulson et al., 2004). 

Jones (2004) found that the average metallicity of planet- 
host stars increased from ~ 0.07 to ~ 0.24 dex for planets 
with semimajor axes of ~ 2 AU to <~ 0.03 AU, suggest- 
ing a trend toward shortest-period planets orbiting the most 
metal-rich stars. Similarly, Sozzetti (2004) showed that both 
metal-poor and metal-rich stars have increasing numbers of 
planets as the orbital period increases but only the metal- 
rich stars have an excess of the shortest period planets. This 
could imply that the metallicity correlation is caused by in- 
ward orbital migration, if low-metallicity stars have long- 
period giant planets that seldom migrate inward. 

Lower disk metallicity results in slower Type II inward 
migration (Livio and Pringle, 2003), the likely dominant 
mechanism for planet migration (see the chapter by Pa- 
paloizou et al). This is because with increased metallicity, 
the disk viscosity v increases. In standard viscous accretion 
disk theory (e.g., Ruden and Pollack, 1991) v = ac s H. 
Lower disk metallicity leads to lower disk opacity, lower 
disk temperatures, lower sound speeds, and a thinner disk. 
As v decreases with lowered metallicity, the time scale for 
Type II migration increases. Ruden and Pollack (1991) 
found that viscous disk evolution times increased by a factor 
of about 20 when v decreased by a factor of 10. It remains 
to be seen if this effect is large enough to explain the rest of 
the correlation. If disk instability is operative and if orbital 
migration is the major source of the metallicity correlation, 
then metal-poor stars should have planets on long-period 
orbits. 

Disk instability may be necessary to account for the 
long-period giant planet in the M4 globular cluster (Sig- 
urdsson et al, 2003), where the metallicity is 1/20 to 1/30 
solar metallicity. The absence of short-period Jupiters in 
the 47 Tuc globular cluster (Gilliland et al, 2000) with 1/5 
solar metallicity could be explained by the slow rate of in- 
ward migration due to the low metallicity. Furthermore, if 



47 Tuc initially contained OB stars, photoevaporation of the 
outer disks may have occurred prior to inward orbital mi- 
gration of any giant planets, preventing their evolution into 
short-period planets, though other factors (i.e., crowding) 
can also be important in these clusters. 

The M dwarf GJ 876 is orbited by a pair of gas giants 
(as well as a much smaller mass planet) and other M dwarfs 
have giant planets as well {Butler et al, 2004), though ap- 
parently not as frequently as the G dwarfs. Laughlin et al. 
(2004) found that core accretion was too slow to form gas 
giants around M dwarfs because of the longer orbital peri- 
ods. Disk instability does not a have similar problem for M 
dwarfs, and disk instability predicts that M, L, and T dwarfs 
should have giant planets. 

With disk instability, one Jupiter mass of disk gas has 
at most ~ 6M ffi of elements suitable to form a rock/ice 
core. The preferred models of the Jovian interior imply that 
Jupiter's core mass is less than <~ 3M ffi (Saumon and Guil- 
lot, 2004); Jupiter may even have no core at all. These mod- 
els seem to be consistent with formation by disk instability 
and inconsistent with formation by core accretion, which 
requires a more massive core. As a result, the possibility of 
core erosion has been raised (Saumon and Guillot, 2004). If 
core erosion can occur, core masses may lose much of their 
usefulness as formation constraints. Saturn's core mass ap- 
pears to be larger than that of Jupiter {Saumon and Guillot, 

2004) , perhaps ~ 15M®, in spite of it being the smaller gas 
giant. Core erosion would only make Saturn's initial core 
even larger. Disk instability can explain the larger Satur- 
nian core mass (Boss et al, 2002). Proto-Saturn may have 
started out with a mass larger than that of proto-Jupiter, but 
its excess gas may have been lost by UV photoevaporation, 
a process that could also form Uranus and Neptune. Disk 
instability predicts that inner gas giants should be accom- 
panied by outer ice giant planets in systems which formed 
in OB associations due to strong UV photoevaporation. In 
low-mass star-forming regions, disk instability should pro- 
duce only gas giants, without outer ice giants. 

Disk instability predicts that even the youngest stars 
should show evidence of gas giant planets (Boss, 1998b), 
whereas core accretion requires several Myr or more to 
form gas giants (Inaba et al, 2003). A gas giant planet 
seems to be orbiting at <~ 10 AU around the 1 Myr-old 
star CoKu Tau/4 (Forrest et al, 2004), based on a spec- 
tral energy distribution showing an absence of disk dust in- 
side 10 AU (for an alternative perspective, see Tanaka et al, 

2005) . Several other 1 Myr-old stars show similar evidence 
for rapid formation of gas giant planets. The direct detec- 
tion of a possible gas giant planet around the 1 Myr-old star 
GQ Lup (Neuhauser et al, 2005) similarly requires a rapid 
planet formation process. 

We conclude that there are significant observational ar- 
guments to support the idea that disk instability, or perhaps 
a hybrid theory where core accretion is accelerated by GI's, 
might be required to form some if not all gas giant planets. 
Given the major uncertainties in the theories, observational 
tests will be crucial for determining the relative proportions 
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of giant planets produced by the competing mechanisms. 
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